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ABSTRACT 

We study, through 2D hydrodynamical simulations, the feedback of a starburst on 
the ISM of typical gas rich dwarf galaxies. The main goal is to address the circulation 
of the ISM and metals following the starburst. We assume a single-phase rotating ISM 
in equilibrium in the galactic potential generated by a stellar disk and a spherical dark 
halo. The starburst is assumed to occur in a small volume in the center of the galaxy, 
and it generates a mechanical power of 3.8 x 10 39 erg s" 1 or 3.8 x 10 40 erg s -1 for 
30 Myr. We found, consistently with previous investigations, that the galactic wind is 
not very effective in removing the ISM. The metal rich stellar ejecta, instead, may be 
efficiently expelled from the galaxy and dispersed in the intergalactic medium. 

Moreover, we found that the central region of the galaxy is always replenished 
with cold and dense gas after a few 100 Myr from the starbust, achieving the requisite 
for a new star formation event in rs 0.5 — 1 Gyr. The hydrodynamical evolution of 
galactic winds is thus consistent with the episodic star formation regime suggested by 
many chemical evolution studies. 

We also discuss the X-ray emission of these galaxies and find that the observable 
(emission averaged) abundance of the hot gas underestimates the real one if thermal 
conduction is effective. This could explain the very low hot gas metallicities estimated 
in starburst galaxies. 
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. 1 INTRODUCTION 



X 



Many gas rich dwarf galaxies are known to be in a starburst 
phase, or are believed to have experienced periods of intense 
star formation in the past (e.g. Gallagher & Hunter 1984; 
Thuan 1991; Tosi 1998, and references therein). These galax- 
ies, classified as "blue compact dwarf (BCD) galaxies" or 
"HII galaxies" , are thus excellent laboratories to investigate 
the feedback of vigorous star formation on the interstellar 
medium (ISM). 

Massive stars inject enormous amount of energy in the 
ISM through stellar winds and when they explode as type 
II supernovae (SNe); the impact of such an energy input on 
the galactic ISM may, in principle, be devastating. In fact it 
is often found that the total energy released during a star- 
burst is greater that the gas binding energy. Yet many dwarf 
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galaxies in a post-starburst phase are still gas rich. As Skill- 
man & Bender (1995) pointed out, observational evidence 
(e.g. Marlowe et al. 1995; Martin 1996) are still insufficient 
to substantiate a disruptive impact of galactic winds on the 
ISM. Clearly, simple energetic considerations do not catch 
the essential nature of the feedback process, and detailed, 
time dependent hydrodynamical models are needed. 

Galactic winds are thought to have a key role in the for- 
mation and evolution of dwarf galaxies (Dekel & Silk 1986, 
Babul & Rees 1992, Matteucci & Chiosi 1983). In general, 
understanding the physics of the feedback of massive stars 
on the ISM is a key problem in cosmological theories of 
galaxy formation (Yepes et al. 1997, Cole et al. 1994). Gas 
outflows from dwarf galaxies are also suggested to be an im- 
portant factor for the production and enrichment of the in- 
tergalactic medium (Trentham 1994). However, persuasive 
arguments against this conclusion are given by Gibson & 
Matteucci (1997), and the origin of metals in clusters of 
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galaxies is still a matter of debate (Brighenti & Mathews 
1998). 

The fate of the (metal rich) material ejected by massive 
stars is of crucial importance in understanting the chem- 
ical evolution of these galaxies (Tosi 1998), in particular 
the low a-elements abundance and the 'strange' values of 
(He/H) and (N/O) vs. (O/H). These problems have been 
encompassed invoking a 'differential ejection', in which the 
enriched gas lost by massive stars escapes from the galaxy 
as galactic wind, while some (or most) of the original ISM 
is unaffected. 

Recent hydrodynamical simulations have verified that, 
under many circumstances, galactic winds are able to eject 
most of the metal rich gas, preserving a significant fraction of 
the original ISM (MacLow & Ferrara 1998, hereafter MF; De 
Young & Heckman 1994; De Young & Gallagher 1991). Silich 
& Tenorio-Tagle (1998) and Tenorio-Tagle (1996), instead, 
found that even the metal-rich material is hardly lost from 
galaxies, since it is at first trapped in the extended halos 
and then accreted back onto the galaxy. 

To investigate this subject further, we present here new 
high resolution calculations, addressing the ultimate fate of 
the ISM and SN ejecta, and their mixing, in a realistic star- 
bursting dwarf galaxy. We investigate in detail the differ- 
ent phases of the gas flow, with particular emphasis on the 
late evolution, evolving the simulations for 500 Myr after 
the starburst event. We consider the effect of the dark mat- 
ter, gas rotation, thermal conduction and different starburst 
strengths. We also discuss the X-ray emission and its diag- 
nostic for the abundance of the hot gas, a particularly ex- 
citing topic in view of the forthcoming launch of AXAF and 
XMM. 

We aim at investigating the evolution of galactic winds 
in a general way, without focusing on any specific object. 
Thus, we select the parameters of the galactic models (total 
mass, ISM mass and distribution, etc.) to be representative 
of the class of dwarf galaxies. Nevertheless, it can be use- 
ful to compare some of our results to a real, representative 
object. An ideal galaxy is NGC 1569, a nearby, well studied 
starburst galaxy. 

Several independent lines of evidence indicate that NGC 
1569 is in a post-starburst phase (Israel 1988, Israel & de 
Bruyn 1988, Waller 1991, Heckman et al. 1995, Greggio et 
al. 1998), with the major starburst activity ceased ~ 5 — 10 
Myr ago. Ha observations of NGC 1569 show (young) bub- 
bles complexes, filaments and arcs throughout the volume 
of the galaxy (Tomita, Ohta & Saito 1994) , suggesting a dif- 
fuse star formation. Heckman et al. (1995) found that the Ha 
emission of NGC 1569 can be separated in a quiescent com- 
ponent, permeating the starbursting region of the galaxy, 
and a more violent component, far more extended and with 
velocities up to 200 km s _1 . This high velocity component is 
interpreted to be ionized shells of superbubbles and provides 
a direct evidence of a galactic-scale outflow. 

Heckman et al. (1995) and Delia Ceca et al. (1996) de- 
tected X-ray emission, extending for 1-2 kpc along the opti- 
cal minor axis of NGC 1569, thus probing the hot gas phase 
directly. This hot gas (T ~ 10 7 K) is the signature of the 
violent SN activity on the ISM. 

As in almost all studies to date, we make a number 
of simplifying assumptions in calculating our models. First, 
the ISM is assumed to be homogeneous and single phase. 



Second, we neglect the selfgravity of the gas, even if the 
gas mass is of the same order of the stellar mass. Third, the 
starburst is instantaneous and concentrated in a small region 
at the center of the galaxy. While none of these hypotheses 
is likely to be strictly correct, they allow for a more direct 
comparison with previous works, and still make possible the 
calculation of models retaining the basics attributes of real 
galactic winds. We will relax some of these assumptions in 
a future paper in preparation. 



2 GALAXY MODELS 

Many ingredients play an important role in determining the 
hydrodynamical evolution of the galactic wind. Among oth- 
ers, the density distribution of the ISM in the pre-burst 
galaxy, the energy injection rate of the newly formed stars, 
the gravitational potential of the galaxy and the effective- 
ness of transport processes in the gas, like thermal conduc- 
tion. 

A thorough exploration of the parameter space would 
require an enormous amount of computational resources and 
it is beyond the scope of this paper. Thus, we hold approx- 
imately constant the stellar and ISM masses of the model 
galaxies (M» = 1.7 x 10 s M Q and M 1S m ~ 1.3 x 10 8 M ), 
although Mism is a crucial factor for the late evolution of 
the system (De Young & Heckman 1994; MF). Instead, we 
vary some of the others parameters as described below. 



2.1 The gravitational potential and the gas 
distribution 

The gravitational potential for our standard model is due 
to two mass distributions: a spherical quasi-isothermal dark 
matter halo plus a stellar thin disk. 

The halo density is given by p h (r) = p h/[l + {r/r c ) 2 ], 
and we chose a central density poh = 4.34 x 10 - 25 g cm -3 
(6.4 x 10 -3 M Q pc" 3 ). The halo core radius is assumed to 
be r c — 1 kpc. The dark halo is truncated at r — 20 kpc. 
The total dark matter mass is thus ~ 2 x 10 9 Mq , while the 
halo mass inside the galactic region (defined hereafter as a 
cylinder R < 2.2 kpc and \z\ < 1.1 kpc, approximately the 
optical size of NGC 1569) is only 0.66 x 10 8 Mq . 

For simplicity, we assume that the stars are distribuited 
in an infinitesimally thin Kuzmin's disk with surface density 



E.(H) = 



2^(^+^)3/2 



where r* = 2 kpc is the radial scalelength and M* = 1.7 x 
10 8 Mq is the total stellar mass, a typical value for dwarf 
galaxies. Although this mass distribution is clearly a rough 
approximation of real stellar disks, it does not degrade the 
accuracy of the large scale hydrodynamical flow. The stellar 
potential generated by this mass distribution is 



$>*{R,z) 



GM, 



^/R 2 + (r* + \z\ 



(Binney & Tremaine 1987). 

It turns out that the stellar mass inside the galactic 
region is M» iga j ~ 3.13 x 10 7 Mq , about half of the dark 
halo mass and about a factor of four less than the gas mass 
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inside the same volume (see below). The dark halo totally 
dominates the mass budget at larger radii. 

The ISM is assumed to be single-phase and in equilib- 
rium with the potential described above. In real dwarf galax- 
ies the neutral ISM is supported against gravity partly by ro- 
tation and partly by the HI velocity dispersion (see Hoffman 
et al. 1996), with maximum rotational velocity that typically 
exceeds the velocity dispersion by a factor of few. Thus, in 
the standard model (hereafter model STD) we allow the 
ISM to rotate, to investigate the role played by the angular 
momentum conservation on the late phase of the evolution, 
when (once the energy output is ceased) the gas tends to 
recollapse toward the central regions (see section 3.1). The 
temperature of the unperturbed ISM is set to To = 4.5 x 10 3 . 

To build a rotating ISM configuration in equilibrium 
with the given potential, we first arbitrarily assume a gas 
distribution in the equatorial plane (z = 0) of the form 
p{R,0) = po/[l + (R/R c ) 2 ] 3/2 , where the central value is 
po = 3.9 x 10 -24 g cm -3 and the gas core radius is R c = 0.8 
kpc. The rotational velocity in the equatorial plane is then 
determined from the condition of equilibrium: 



2 

v <t> 



dp 



dR 



where v c = y Rd^/dR is the circular velocity and p the 
thermal gas pressure. The rotational velocity is assumed to 
be independent of z. The density at any z is then found 
integrating the z-component of the hydrostatic equilibrium 
equation, for any R. The edge-on and face-on profiles of the 
resulting gas column density are shown in Fig. 1. We note 
that this model, having an extended gaseous halo, resembles 
the models worked out by Silich & Tenorio-Tagle (1998). 

The circular velocity for this mass model increases with 
R, reaching the maximum value of ~ 20 km s _1 at R ~ 4 
kpc and staying almost constant for larger 7?. The rota- 
tional velocity v<f, shows a similar radial behaviour, but with 
a maximum value ^ 15 km s . 

The total gas mass inside the galactic region is 
A^iSM, g ai ~ 1-32 x 10 s Mq , a typical amount for dwarf 
galaxies (Hoffman et al. 1996), and in close agreement with 
the mass inferred for NGC 1569 in particular (Israel 1988). 
The total gas mass present in the numerical grid (extending 
to 25 kpc in both R and z directions) is Misin.tot ~ 6 x 10 8 
M . 

The gas distribution qualitatively resembles that used 
by Tomisaka & Ikeuchi (1988). It has a low density region 
around the z-axis (see Fig. 2a), which acts as a collimat- 
ing funnel for the hot outflowing gas (Tomisaka & Bregman 
1993, Suchkov et al. 1994). This is due to the assumption 
that Vj, does not depend on z. The funnel, however, influ- 
ence the gas dynamics only at very large distances above the 
galactic plane (i.e. for z > 10 kpc) and does not invalidate 
the results presented in sections 3 and 4. 

In order to address the influence of an intracluster 
medium (ICM) confining the galactic ISM, we calculate 
model PEXT (section 4.3). In this simulation we replace 
all the cold ISM (distributed as described above) having 
a thermal pressure P < 10~ 13 dyn cm -2 with a hot, rar- 
efied ICM with picm = 8 x 10" 30 g cm -3 and Ticm = 10 8 
K. In this case the cold ISM is confined to a roughly el- 
lipsoidal region with major and minor semiaxes 2 kpc and 
1 kpc, respectively. The galactic ISM mass is now only 



Figure 1. Column density of the initial ISM. Heavy solid line: 
model STD seen edge-on; heavy dashed line: model STD seen 
face-on; light solid line: model B seen edge-on; light dashed line: 
model B seen face-on. 



A^iSM, g ai = 1-05 x 10 s Mq , while the total mass of gas 
in the grid is MisM.tot ~ 1-16 x 10 8 M Q . 

In addition to the models described above, we use a dif- 
ferent galaxy model (model B) to investigate the effect of the 
absence of dark matter and rotation. The isothermal ISM in 
hydrostatic equilibrium in the potential well generated by 
the same stellar distribution as in model STD. The central 
gas density is po = 1.1 x 10~ 23 g cm -3 , and the gas mass 
inside the galactic region is 1.4 x 10 8 Mq , approximately 
as in model STD. Due to the lack of rotational support, the 
gas distribution is now more concentrated than in model 
STD (see the ISM column density in Fig. 1), and the total 
gas mass inside the grid is MrsM.tot = 2.3 x 10 8 Mq . We 
also run a model identical to model B, but including heat 
conduction (model BCOND). 



2.2 The starburst 

We assume an instantaneous burst of star formation which 
injects energy in the ISM for a period of 30 Myr, approx- 
imately the lifetime of a 8 Mq star, the smallest star pro- 
ducing a type II SN. We consider two starburst strengths: 
the first is representative of a moderate starburst, while the 
second is intended to match more active galaxies. We as- 
sume that the starburst produces a steady (mechanical) en- 
ergy input rate L- mp = 3.76 x 10 39 erg s _1 (hereafter SB1 
model) and L inp = 3.76 x 10 40 erg s" 1 (SB2 model) respec- 
tively. Model SB2 produces a mechanical power similar to 
the lower limit estimated for NGC 1569 (Heckman et al. 
1995). The mass injection rate is assumed to be respectively 
M = 3 x 10~ 3 M yr" 1 and M = 3 X 10~ 2 Mq yr" 1 . 

It is useful to compare the parameters we use with the 
detailed starburst models of Leitherer & Heckman (1995) 
(LH). For example, for an instantaneous burst with a 
Salpeter IMF (from 1 to 100 M ) and metallicity 1/4 of 
solar, they found a mechanical energy deposition rate ap- 
proximately constant between 6 and 30 Myr after the star- 
burst (see their fig. 55). This justifies our assumption of a 
steady energy source. Our assumed mechanical luminosities 
for SB1 and SB2 correspond respectively to ~ 2.1 x 10 and 
~ 2.1 x 10 6 Mq turned into stars during the starburst event, 
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Table 1. Models physical and numerical parameters. 
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405 2 


480x540 


410x530 


A/? ■ 


10 


3 


10 


2 


2 


A 7?„, 


zo 


1 7 






99 


^max 


25 


14.9 


25 


23 


11.4 


2max 


25 


14.9 


25 


42 


38 


Code 


ZEUS 


ZEUS 


ZEUS 


BO 


BO 



A^ISM is the initial ISM mass present in the computational grid; 
MisM,gal is the ISM mass in the galactic region R < 2.2 kpc and 
\z | < 1.1 kpc; M, is the assumed stellar mass; Mdh is the dark 
matter halo mass. All masses are given in units of 10 8 Mq . Lj np is 
the energy input rate in 10 39 erg s — 1 . and N z are the number 
of cells in the ij-direction and z-direction. A/? m ; n (=Az m j n ) is 
the central zone size in pc. Ai?2k pc (=^ 2 2kpc) is the width in pc 
of the zone at (R,z)=(2 kpc,0) (or (R,z)=(0,2 kpc)). i? max and 
z m ax are the total dimensions of the grid in kpc. In the last row 
is indicated the hydrocode used (ZEUS-2D or the Bologna code). 

according to LH. If all stars with initial mass greater than 8 
Mq end their lives as type II supernovae, the total number 
of SNII events produced by the starbursts is ~ 4000 and 
~ 40000 for SB1 and SB2. 

The total energy deposited after 30 Myr is ~ 3.56 x 10 54 
erg and ~ 3.56 x 10 55 for SB1 and SB2. These values must be 
compared with the binding energy of the gas present in the 
numerical grid in the standard model, -Ebmd ~ 5.3 x 10 54 erg, 
and with the binding energy of the gas inside the galactic 
region ~ 1.8 x 10 54 erg. 

After 30 Myr, the total mass returned to the ISM by 
stellar winds and type II SNe in the model by LH, again 
with Z =1/4 Z Q and a Salpeter IMF, is 4.81 x 10 4 Mq and 
4.81 x 10 5 M for SB1 and SB2. With our assumed M, 
however, we inject 9 x 10 4 M and 9 x 10 5 Mq for models 
SB1 and SB2 respectively. Thus, we overestimate the mass 
return rate by a factor ~ 2 with respect to the LH model ^} 
However, the hydrodynamical evolution of our models is not 
sensitive to such a discrepancy, as well as our estimate of the 
efficiency of ISM and metal ejection (although the pollution 
degree of the ISM may be affected). 

While our assumed starburst model is fairly consistent 
with the detailed theoretical models by LH, it is impor- 
tant to note that real galaxies have generally a much more 
complex star formation history. The assumption of instanta- 
neous, point-like burst appears particularly severe. For ex- 
ample, Greggio et al (1998) found that the bulk of the star- 
burst in NGC 1569 proceeded at an approximately constant 
star formation rate of 0.5 Mq yr _1 for 0.1 — 0.15 Gyr (as- 
suming a Salpeter IMF from 0.1 to 120 Mq ), until ~ 5 - 10 
Myr ago, when the star formation in the field ended. It im- 
plies that ~ 5 — 7.5 x 10 7 Mq of gas has been converted into 
stars. 

t Alternatively, we can think to a starburst with a double amount 
of mass turned into stars, and to an efficiency in the energy de- 
position rate of ~ 0.5. 



Moreover, Hcv observations of NGC 1569 show 
(young) bubbles complexes, filaments and arcs distributed 
throughtout the volume of the galaxy (Tomita et al. 1994), 
suggesting a diffuse, wide scale star formation. 

Galactic wind models powered by a point-like energy 
source are nevertheless useful as first step toward the full 
complexity of the problem, and for a direct comparison with 
previous studies. Simulations with spatially and temporally 
extended star formation will be the subject of a forthcoming 
paper. 

2.3 The numerical simulations 

To work out the models presented in this paper we used two 
different 2-D hydrocodes. The first one has been developed 
by the Numerical Group at Bologna Astronomical Observa- 
tory and the (1-D) core of the scheme is described in Bedogni 
& D'Ercole (1986). This code and its successive extensions 
have been applied to a variety of astrophysical problems 
(e.g. Brighenti & D'Ercole 1997, D'Ercole & Ciotti 1998). 
The second code employed is ZEUS-2D, a widely used, well 
tested scheme developed by M. Norman and collaborators 
at LCSA (Stone & Norman 1992). We always found consis- 
tent results among the codes, as expected from the numer- 
ous hydrodynamical tests performed with the Bologna code 
(Brighenti 1992). 

We solve the usual hydrodynamical equations, with 
the addition of a mass source term and a thermal energy 
source term; the hot gas injected expands to form the star- 
burst wind with the appropriate mechanical luminosity L[ np . 
These equations are described in details in, e.g., Brighenti 
& D'Ercole (1994). The (constant) mass and energy source 
terms are given respectively by a — M /V and at. Here V 
is the volume of the source region, chosen to be a sphere of 
radius=50 pc, centered at (R,z) = (0,0), and e = L lnp /M. 
To keep track of the gas lost by the stars formed in the 
starburst (the ejecta), we passively advect it solving an an- 
cillary continuity equation for the ejecta density p c j. Both 
the codes used spread shocks over 3-4 zones and contact 
discontinuities over 4-10 zones. 

In our models the angular momentum is treated in a 
fully consistent way (see Stone & Norman 1992 for the 
details about the resolution of the angular momentum 
equation). Thus, contrary to some of the previous studies 
(Tomisaka & Ikeucki 1988, Tomisaka & Bregman 1993), we 
do not use a reduced gravitational force to mimic the rota- 
tional support of the ISM. 

To take into account the thermal conduction (model 
BCOND) we adopt the operator splitting method. We iso- 
late the heat diffusion term in the energy equation and solve 
the heat transport equation, alternatively along the z and R 
direction separately, through the Crank-Nicholson method 
which is unconditionally stable and second order accurate. 
The system of implicit finite difference equations is solved 
according to the two-stage recursion procedure (e.g. Rich- 
meyer and Morton 1967). Following Cowie & McKee (1977), 
we adopt saturated fluxes to avoid unphysical heat transport 
in presence of steep temperature gradients. 

We ran the models on a cylindrical grid (coordinates R 
and z), assuming axial symmetry. We use reflecting bound- 
ary conditions along the axes and outflow boundary condi- 
tions at the grid edges. To better resolve the central region, 
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the grid is unevenly spaced, with the zone width increas- 
ing from the center to large radii. Specifically, in the stan- 
dard model (STD+SB2), the grid extends in both the R- 
direction and z-direction from to 25 kpc. The first zone is 
AT? = Az = fO pc wide, and the size ratio between adjacent 
zones is 1.00747. 

For the other models we use different grid spacing. For 
model SB1 (§4.1) the inner grid size is 3 pc and the size 
ratio is 1.00717. For model B and BCOND the central zone 
is only 2 pc wide and the size ratio between adiacent zones 
is 1.01. 

The parameters used in the models and other charac- 
teristic quantities are summarized in Table 1. 



3 THE STANDARD MODEL (STD+SB2) 
3.1 The dynamics of the ISM 

As the starburst wind starts blowing, the classical two 
shocks configuration is achieved, in perfect analogy to the 
stellar wind bubble theory (Dyson & de Vries 1972, Weaver 
et al. 1977). The freely expanding wind encounters the re- 
verse shock and is heated to T ~ 5 x 10 7 K, while the ex- 
ternal shock sweeps the ISM. The shocked starburst wind 
and the shocked ISM are separated by a contact disconti- 
nuity. The reverse shock is always approximately spherical, 
since the short sound crossing time in the shocked wind re- 
gion keeps the pressure almost uniform. The shape of the 
forward shock, instead, depends on the ISM density distri- 
bution. The roughly oblate ISM configuration of our model 
galaxy forces the superbubble to expand faster along the 
polar direction, acquiring the classical bipolar shape at late 
times (Fig. 2b,c,d). At earlier times, however, the density 
distribution favours a diagonal expansion, generating a cu- 
rious boxy morphology (Fig. 2a). 

When the energy input from the starburst ends (at 
t = 30 Myr, Fig. 2a), almost the whole galactic region 
(R < 2.2 kpc, \z\ < 1.1 kpc) is filled by the freely expanding 
wind, a situation clearly unrealistic, due to our simple ISM 
model. In real galaxies we expect that this region hosts a 
complicated multiphase medium. Israel & van Driel (1990) 
found a relatively small hole (with diameter ~ 200 pc) in 
the HI distribution, associated with the 'super star cluster' 
N 1569A and likely caused by the action of SNII and stellar 
winds of the star cluster. 

The shocked ISM shell is accelerating through the steep 
density gradient of the unperturbed ISM, and this accelera- 
tion promotes Rayleigh- Taylor (R-T) instabilities. The shell 
tends to fragment and relatively dense (n < 0.5 cm -3 ), 
cold (T ~ 10 4 K) filaments are clearly seen in Fig 2a at 
(z,R) ~ (1,3) kpc. The numerical resolution of this simula- 
tion is not appropriate to follow the real formation and evo- 
lution of these features, whose actual density is expected to 
be higher than that found in our computations. As pointed 
out by MF, these filaments of dense gas are immersed in the 
hot gas and they do not necessarily trace the edges of the 
superbubble. 

At t=30 Myr the whole shocked ISM shell is radiative. 
Following the trend of the external ISM, a density gradient 
is present along the shell, from the equator to the pole, with 
the lateral side being denser (n ~ 1.5 cm" 3 ) and the polar 



Figure 2. Map of the logarithm of the number density (cm -3 ) 
for model STD at four different times. The z-axis is vertical. Dis- 
tances are in kpc. The gray-scale varies linearly from -7.65 (white) 
to 0.135 (black). Six contours are shown at levels -1, -2, -3, -4, -5 
and -6. 



region more rarefied (n ~ 0.002 cm -3 ). At this time, the 
optical appearance of the system, if kept ionised (for exam- 
ple by a low level star formation activity, hot evolved stars, 
etc.), would be that of an incomplete shell, the polar portion 
being too rarefied to be observable, given its low emission 
measure (EM « 0.01 cm" 6 pc for the shell and EM < 30 
cm -6 pc for the filaments). 

The radiative external shock is too slow (v s hock ^ 180 
km s -1 ) to emit X-ray. This is contrary to the results by 
Suchkov et al. (1994), who claim that most of the X-ray 
radiation comes from the shocked ISM. The different be- 
haviour in our models is due to the lower (by an order of 
magnitude) Li np considered, probably more appropriate for 
a dwarf galaxy. As explained in section 3.4, we also find that 
the hot ISM is the most important contributor to the X-ray 
luminosity. However, it is not heated by shocks, but by mix- 
ing with the shocked wind at the contact discontinuities. At 
later times the outer shock accelerates through the steep 
density gradient, and heats the ISM to X-ray temperatures. 
On the other hand, this happens only when the X-ray lumi- 
nosity has dropped to very low and uninteresting values (cf. 
Fig. 6). 

Fig. 2b shows the density distribution at t = 60 Myr. 
The steep density gradient along the ^-direction induces a 
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radiative-adiabatic transition of the polar portion of the ex- 
ternal shock. The cold filaments are slowly moving forward, 
and their density decreases to maintain the pressure equilib- 
rium with the expanding hot gas; now the densest filaments 
have n ~ 0.04 cm -3 . 

The shell is increasingly thicker with time. In fact, while 
the outer edge is still expanding, with v ~ 25 km s , the 
inner edge of the shell near the equatorial plane is already 
receding toward the center with a velocity v ~ 10 km s _1 . 
This backward motion, due to the drop of the pressure inside 
the expanding hot bubble, will eventually cause the collapse 
of the cold gas back inside the galactic region, as evident in 
Fig. 2d (see also MF) . The details of the collapse are shown 
in Fig. 3 and described below. 

Fig. 2c and 2d show the density at 100 and 200 Myr re- 
spectively. The external shock assumes a prounounced cylin- 
drical shape because of the collimating effect of the low den- 
sity region around the z-axis (section 2.1). The polar shock 
crossed the numerical grid edge (at z — 25 kpc); however, 
being the motion supersonic, the numerical noise generated 
at the grid boundary does not propagate back. Moreover, 
given the very low densities in that region, the amount of 
gas lost from the grid is completely neglibile. 

The temperature of the hot, X-ray emitting gas de- 
creases with time. T ~ 5 x 10 7 K during the active energy 
injection phase (t < 30 Myr). At later times, radiative losses 
and especially expansion lower the temperature (T < 10 7 K 
at 60 Myr; T < 2 x 10 6 K at 100 Myr; T < 10 6 K at 200 
Myr). ASCA observations of NGC 1569 (Delia Ceca et al 
1996) indicate that the diffuse X-ray emission comes from a 
(luminosity weighted) T ~ 8 x 10 6 K gas. We note, however, 
that the temperature inferred from simple fits to observed 
X-ray spectra may be a poor estimate of the actual tem- 
perature (Strickland & Stevens 1998). We warn that all the 
model temperature values quoted are mass weighted and 
may not represent the "observable" ones. 

No more supported by the hot gas pressure, the cold 
ISM recollapses at t ~ 150 Myr, filling again the galactic 
region. We show a zoomed view of the central part of the grid 
in Fig. 3. In panel a the density contours and the velocity 
field are shown at t = 140 Myr, just before the cold gas 
reaches the center. In panel b the same quantities are shown 
at t — 200 Myr. In Fig. 3a the cold tongue is approaching 
the center at v ~ 40 km s _1 , with a Mach number M ~ 10. 
The density in the collapsing gas increases with R, from n ~ 
2.3 x 10" 3 cm" 3 at (R,z) = (0.5,0) kpc, to n ~ 7.5 x 10" 3 
cm" 3 at (R,z) = (1,0) kpc, ton ~ 2.3 x 10~ 2 cm" 3 at 
(R, z) = (2, 0) kpc. At t = 150 Myr the cold gas reaches the 
center and shocks. Hereafter the accretion proceeds through 
a cylindrical shock wave. 

At t = 200 Myr the cold gas is still flowing toward the 
center with v ~ 15 km s _1 , building a (transient) conical 
structure around the z-axis (Fig. 3b). The mean density on 
the equatorial plane in the galactic region is n ~ 0.025 cm -3 , 
and it is still growing with time. Panels c and d show the 
subsequent evolution (t = 250 Myr and t = 300 Myr respec- 
tively). The collapse is the result of the pressure drop in the 
hot bubble, due to its expansion along the polar direction. 
Between t = 60 Myr and t = 100 Myr the pressure of the hot 
gas decreases by more than one order of magnitude (from 
~6x 10" 13 to ~ 10" 15 dyn cm" 2 ). The cold gas, no longer 



supported by the hot phase, is driven back to the center 
mainly by its own pressure, rather than the galactic gravity. 

We have verified that the collapse is not a spurious nu- 
merical effect due to strong radiative losses at the numeri- 
cally broadened contact discontinuities. To this purpose we 
run a simulation identical to the standard model, but with 
the radiative cooling turned off. For this adiabatic model 
we found that the collapse time is again ~ 150 Myr. The 
ISM recollapse is thus an unavoidable phenomenon for all 
the models considered in this paper. 

It is interesting to follow the circulation of the gas. The 
mass of the ISM in the galactic region, at t = 200 Myr, 
is MisM.gai ~ 0.034 x 10 s M , about a factor of ~ 40 
lower than the initial gas mass in the same volume (note 
that it does not mean that the ejection efficiency /ism is 
1/40, since what matters in estimating /ism is the amount 
of ISM effectively bound to the galaxy, as described in the 
next section). However, at t = 200 Myr the cold gas is still 
flowing toward the center (Fig. 3b); for instance, at t = 300 
Myr we found MisM, ga i = 0.096 x 10 s Mq . At t = 500 Myr, 
the final time of our simulation, MisM.gai = 0.27 x 10 s Mq , 
a factor of ~ 5 lower than the initial gas mass in the same 
volume. 

We can speculate further on the fate of the cold gas 
falling back to the center at late times. The face-on sur- 
face density of the central ISM is an increasing function of 
time (Fig. 4), since material continues to accrete until the 
end of our simulation (t = 500 Myr; we could not follow the 
evolution further because of our limitation in computational 
resources). It has been suggested that above a critical ISM 
surface density £ C rit ~ 5 — 10 x 10 20 cm" 2 , the star forma- 
tion in dwarf galaxies is very efficient (Gallagher & Hunter 
1984, Skillman 1987, 1996). At t = 200 Myr the face-on sur- 
face density peak in our model is E ~ 10 20 cm -2 , and it 
grows slowly to 6 x 10 20 cm -2 at t = 500 Myr. Thus, we can 
hypothesize that the threshold surface density is reached in 
a time of the order of 1 Gyr, after that a new burst of star 
formation may start. This scenario is rougly consistent with 
many studies of the star formation history in BCD galaxies, 
which indicate that stars are formed mainly through several 
discrete, short bursts, separated by long (~ few Gyrs) qui- 
escent periods (see the review by Tosi (1998) and references 
therein). 

3.2 The ISM ejection efficiency 

A key point in the galactic wind theory is the ability of the 
starburst in ejecting the ISM (see Skillman & Bender (1995) 
and Skillman (1997) for a critical review about this subject). 

We estimate the ISM ejection efficiency calculating, at 
some late time, for example t=200 Myr, the mass of ISM 
Miost which has velocity or sound speed greater than the 
local escape velocity. We assume that this gas (and the gas 
that already left the grid) will be lost by the galactic system 
(see also MF). It is important to note that Mi ost calculated 
in this way should be considered only a rough estimate of 
the amount of gas leaving the galaxy, since dissipative effects 
may lower the ejection efficiency, and the escape velocity de- 
pends critically on the poorly known size of real dark matter 
halos. The ISM ejection efficiency, /ism, is then defined as 
Miost /A/initial) where iWinitiai is the total gas mass present 
on the whole grid at t — (we neglect the contribution of 
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Figure 3. Density contours and velocity field for the central region of model STD at 140 Myr, 200 Myr, 250 Myr and 300 Myr. 



thus, the efficiency in removing the ISM from the central 
regions may be considerably greater than /ism- However, for 
other models (see section 4.3), the galaxy is able to recover 
most of the original ISM mass after pb 100 Myr. 



Figure 4. Face-on column density of model STD at three differ- 
ent times. Dotted line: 200 Myr; dashed line: 300 Myr; solid line: 
500 Myr. 

the ejecta, whose total mass is only < 0.2 % of the initial 
mass) . We note that this operative definition for /ism is grid 
dependent, since Mi n i t i a i increases with the volume covered 
by the numerical grid. 

At t=200 Myr we find /ism = 0.058: evidently even a 
powerful starburst as the one considered for this model is 
not effective in removing the interstellar gas. 

However, as pointed out in the previous section, the 
gas mass inside the galactic region can be significantly lower 
than the initial value, even long after the starburst event: 



3.3 The enrichment 

In a similar way we have estimated the ejection efficiency of 
the metal-rich stellar ejecta, / c j. At t = 200 Myr we found 
/ c j = 0.46: the galaxy is less able to retain the enriched stel- 
lar ejecta than its own original ISM. This finding supports 
the selective winds hypothesis, and it is in qualitative agree- 
ment with others numerical simulations (De Young & Gal- 
lagher 1990; De Young & Heckman 1994; MF). 

It is interesting to investigate the spatial distribution 
of the ejecta material. We found that at t = 200 Myr 
~ 7.3 x 10 5 Mq of stellar ejecta are present on the nu- 
merical grid, about 80 % of the total material released by 
the starburst. However, the ejecta mass in the galactic re- 
gion is only M ojiga i ~ 5.15 x 10 3 M Q , less than 0.6 % of 
the total amount ejected (9 x 10 5 Mq )! Since gas continues 
to flow toward the central region, the mass of the ejecta in 
the galactic region increases slightly with time. At t = 300 
Myr, for instance, we found M j, ga i ~ 1.29 x 10 4 Mq , and 
M ejiga i ~ 3.6 x 10 4 M Q at t = 500 Myr. We conclude that, 
while a significant fraction of the ejecta is retained by the 
relatively deep potential of the dark halo, most of it resides 
in the outer regions of the system, in a phase so rarefied to 
be virtually unobservable. 

The cold gas collapsing at late times, and filling the 
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galactic region, has been only slightly polluted by the stellar 
ejecta. To characterize the pollution degree we introduce the 
local ejecta fraction as Z = p c j / p, where p c j is the density of 
the ejecta. The average ejecta fraction in the galactic region, 
at t — 200 Myr, defined as < iJ ga i >= M e j !ga i/MisM, g ai, is 
Z ~ 1.4 X 10~ 3 . The cold galactic ISM, probably the only 
component detectable at late times because of its relatively 
high density, shows only a small degree of enrichment. 

We can estimate the increase in the metal abundance 
generated by the starburst from the total number of SNII, 
which we assume to be the only source of metals. The iron 
production and circulation is particularly worthwhile, be- 
cause the metallicity estimated through X-ray spectra of the 
hot gas phase (T ~ 10 7 K) are especially sensitive to iron 
through the Fe-L complex at ~ 1 keV. For the sake of sim- 
plicity we shall neglect the iron produced by SNIa, whose 
iron release timescale is believed to be of the order of one 
Gyr (Matteucci & Greggio 1986), a time much longer than 
those considered in this paper. 

In section 2.2 we estimated a total number of SNII ~ 
4000 and ~ 40000 for SB1 and SB2 respectively (adopting 
the same IMF as in LH). The yields of metals from SNII 
are rather uncertain, especially for iron and oxygen, because 
of the complications in the late evolution of massive stars 
and nuclear reactions rates. A compilation of IMF averaged 
yields, i.e., the mean ejected mass of a given element per 
SN, can be found in Loewenstein and Mushotzky (1996). 
Given the approximate nature of the calculations presented 
in this paper, we simply adopt < j/f c >= 0.1 Mq and < 
yo >= 1 Mq as reasonable values for averaged iron and 
oxygen yields. We assume that the metals are well mixed 
within the ejecta, whose abundances (by mass and relative 
to H) are Zf c ^ ~ 3.4 Zf b ,q and Zo,ej ~ 4.6 Zo,q, where 
we adopt the meteoritic solar abundances from Anders & 
Grevesse (1989). 

In Fig. 5 we show the iron gas abundance distribution 
at t = 200 Myr, assuming that the original ISM has Zf c = 
(i.e. we calculate the increment in the metallicity caused by 
the starburst ejecta). The iron abundance is highly inho- 
mogeneous, both in the hot and cold phase. It ranges from 
very low values Zfc < 0.01 Zf c ,@, to the pure ejecta value 
Z-p c = 3.4 Zf B] q. The hot phase metallicity is supersolar 
with typical values Zp a = 1.5 — 2.5 Zf c ,q- It is puzzling 
that recent ASCA observations of the outflows in starburst 
galaxies indicate that the metal abundance of the hot gas is 
rather low. We discuss this point in the next section. 

While the numerical diffusion may affect somewhat the 
absolute values of Zf c , we believe that the spatial variations 
of the metallicity are real. The oxygen abundance pattern is 
identical to that of iron, due to our assumption of perfectly 
mixed ejecta, but with different minimum and maximum 
values (0 - 4.6 Zo,©). 

The cold gas replenishing the galactic region has average 
metallicity < Z Fc >< 0.005 Z Fe , (< Z >~ 0.01 Zo,©), so 
that a successive instantaneous starburst event would form 
stars only slightly more metallic than the previous stellar 
generation. 

3.4 The X-ray emission and hot gas metallicity 

A detailed investigation of the X-ray emission of the hot gas 
is beyond the scope of this paper. The intrinsic diffusion of 



Figure 5. Distribution of the iron abundance, calculated as ex- 
plained in the text, for model STD at t = 200 Myr. Labels are in 
kpc. Minimum and maximum values, in solar units, are (white) 
and 3.4 (black). Contours are shown at Zp c = 0.5, 1,1.5, 2, 2.5, 3. 



the numerical scheme spreads contact discontinuities, sepa- 
rating the hot and cold phases, over several grid points. This 
fact prevents us to consistently calculate X-ray luminosities 
(Lx) and emission averaged abundances of the hot phase. In 
fact, the gas inside the broadened contact discontinuities (a 
mixture of the ejecta and the pristine ISM), being relatively 
dense and with temperatures of the order of 10 6 K, turns 
out to dominate Lx- 

It is important to note that several physical processes, 
not considered in this simulation, smear out hydrodynam- 
ical discontinuities, mixing cold and hot ISM, and produc- 
ing a gas phase with intermediate temperature and density. 
Most important are thermal conduction and turbulent mix- 
ing (Begelman & Fabian 1990). Thus, the undesired numer- 
ical diffusion qualitatively mimics real physical effects. We 
consider explicitely the heat conduction in section 4.3, but 
we are not in the position to make quantitative estimates on 
the influence of turbulent mixing layers on Lx ■ 

With this limitation in mind we can nevertheless gain 
some insights on the properties of the X-ray emission of 
starbursting galaxies. The X-ray luminosity of several mod- 
els (see below), calculated in the straightforward way as 
Lx = J en(T)dV, is a decreasing function of time (here 
cr(T) is the Raymond-Smith emissivity in the ROSAT 
band). All models shown in Fig. 6 share the same trend: 
at first Lx drops gently until t = 30 Myr and then, when 
the energy input stops, Lx decreases rapidly to unobserv- 
able values. 

One of the most interesting observables is the emission 
averaged metallicity. In order to compare our data with os- 
ervations, we define the emission averaged ejecta fraction 
< Z >x~ (1/Lx) J ZeadV (it is a measure of the gas 
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Figure 6. X-ray luminosity vs. time for four models. Solid line: 
model STD; dashed line: model SB1; dotted line: model B; dot- 
dashed line: model BCOND. 

abundance Q). It is almost constant with time up to 30 
Myr, with typical values of 0.04-0.06, , and then it increases 
steadly up to ~ 0.5 at t = 100 Myr, when, however the 
X-ray luminosity is so weak to make the detection virtually 
impossible. This very low 'metallicity' clearly indicates that 
most of the X-ray emission comes from original ISM mixed 
with stellar ejecta. 

To isolate the contribution of the ejecta material to Lx 
and < Z >x, we recalculated these two quantities using 
the ejecta density p CJ (instead of the gas density) in the 
calculation of the X-ray emissivity. Now Lx is a factor ~ 
100 lower than before, and < Z >x is now much larger, 
approaching unity. Note that this does not mean that the 
ISM heated by the outer shock is the main contributor to Lx 
(see below) . It demonstrates instead that diffusion processes 
mix the ISM with the shocked ejecta, and the material in 
these mixing layers makes most of the Lx ■ 

In the last few years the ROSAT, ASCA and BeppoSAX 
satellites provided detailed X-ray observations of starburst 
galaxies. While for dwarf galaxies the hot gas abundance 
cannot be unambiguosly determined (e.g. Delia Ceca et al. 
1996), for brighter starburst galaxies the counts statistics 
is high enough to make this task possible (Ptak et al 1997; 
Tsuru et al. 1997; Okada, Mitsuda & Dotani 1997; Persic et 
al. 1998) . A somewhat surprising result of all these observa- 
tions is that the iron abundance is invariably small, typically 
less than 0.1 solar. This low metallicity can easily be under- 
stood if the X-ray emission is dominated by the layer of 
shock-heated ISM, as pointed out by Suchkov et al. (1994). 
However, this is not a general result, and it does not hold 
for our models in particular, since the external shock is too 
slow to heat the ISM to X-ray temperatures (section 3.1). 
Thus, in model STD the only X-ray emitting gas is expected 
to be the (shocked) ejecta of the stars formed in the star- 

§ The emission averaged iron abundance is Zp c<c j < Z >x, where 
•^Fe.ej = 3.4 Z-p Ct Q with the assumptions adopted in section 3.3. 



burst, and its metallicity is thus expected to be quite high. 
This abundance discrepancy forces the theoretical models to 
move toward a higher level of complexity. 

Low X-ray abundances can be explained in several ways. 
First, it seems reasonable that thermal conduction and tur- 
bulent mixing give rise to a mass loaded flow (Hartquist, 
Dyson and Williams 1997, Suchkov et al. 1996) with low 
emission averaged metallicity, provided that the cold gas 
mixed with the hot phase is nearly primordial. In this case 
the emission averaged temperature of the hot gas is expected 
to be low (few 10 6 K); see section 4.3. Second, the hot gas 
might be severly depleted by dust. Stellar outflows and SN 
ejecta are observed to form dust (e.g. Clegg 1989; Colgan et 
al. 1994), and the dust sputtering time t Bp ~ 2 x 10 6 a Mm /n 
yr ( where a, lm is the dust grain radius, Draine & Salpeter 
1979; Itoh 1989) in the hot phase may be long enough to 
make most of the iron still locked into grains after few 10 7 yr. 
Another possibility is that the estimated abundances are not 
accurate. Strickland & Stevens (1998) analysed the synthetic 
ROSAT X-ray spectrum of a simulated wind-blown bubble, 
finding that simple fits may underestimate the metallicity by 
more than one order of magnitude. The inadequacy of 1-T 
models in estimating the gas abundance has been demon- 
strated also by Buote & Fabian (1998) and Buote (1999) in 
the context of hot gas in elliptical galaxies and groups of 
galaxies. Indeed, Dahlem, Weaver & Heckman (1998) used 
multi-components models to fit ROSAT PSPC + ASCA 
spectra of seven starburst galaxies and found that low metal- 
licities are no more required, and nearly solar abundances 
are entirely consistent with the data. Their findings support 
the idea that the inferred low abundances are caused by the 
undermodelling of X-ray spectra. 

4 OTHER MODELS 
4.1 Model SB1 

With this simulation we investigate the effect of a weaker 
starburst on the ISM of a dwarf galaxy. This model is iden- 
tical to model STD, but the starburst mechanical luminosity 
is a factor of ten lower. This starburst may be more typical 
among dwarf galaxies. We reduce the mechanical luminosity 
lowering the mass loss rate by a factor of ten (see §2.2). 

We anticipate that in this model the radiative cooling at 
the contact discontinuities, artificially broadened by numer- 
ical diffusion, is now important, and causes the hot bubble 
to slowly collapse. 

Fig 7a shows the density distribution at t = 30 Myr, 
when the energy and mass input turns off. As in model STD, 
dense tongues of shocked ISM penetrate in the hot bubble as 
a result of R-T instabilities. As expected, the superbubble is 
now much smaller than in model STD, and it is expanding 
less rapidly. The hydrodynamical evolution is illustrated in 
the other panels of Fig. 7. At * = 60 Myr (Fig. 7b) the in- 
ternal edge of the shocked ISM shell is receding toward the 
center with v ~ 25 km s _1 . The cold gas reaches the origin 
at t ~ 70 Myr, much earlier than in model STD. We find 
that the ISM collapse is slightly accelerated by the spurious 
energy losses mentioned above. To address the importance 
of this undesired numerical effect, we recalculated the Model 
SB1 without radiative cooling (the adiabatic model). These 
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starburst events in this model is shorter that in model STD, 
being only few 100 Myr. 

At t = 200 Myr the ISM ejection efficiency /ism is es- 
sentially zero: all the gas is cold and it is moving with a 
velocity lower than the escape velocity (/ism = 1.8 X 10~ 3 
for the adiabatic model). The gas mass inside the galactic 
region is MisM.gai = 6.5 x 10 7 Mq , about half of the mass 
present initially. Since the gas is still accreting, the central 
ISM mass increases with time: at t = 300 Myr we have 
MiSM.gai = 8.0 x 10 7 M . 

Thus, in the case of moderate starburst strength, the 
galaxy is able to recover most of the original ISM in a rela- 
tively short time. The evolution of this model is qualitatively 
similar to that of model STD, but is now accelerated and, as 
expected, the galactic ISM 'forgets' the starburst quicker. 

The circulation of the stellar ejecta is qualitatively simi- 
lar to that of the standard model. However, now / e j = 0.003: 
almost all the metals produced by the starburst remain 
bound to the galaxy. A significant fraction of the total ejecta 
mass (~ 2.4 x 10 4 M , ~ 27 % of the total) is still present 
in the galactic region at this late time. The very low value 
for / e j is partly due to the excess of radiative losses at the 
contact surfaces. For the adiabatic model we find / e j = 0.14 
(still much lower than in model STD) and M c j iga i ~ 3.3 x 10 4 
Mq . In summary, we find that / e j is significantly lowered 
by the spurious extra-cooling, but the important quantity 
Afej.gai does not change greatly. The conclusion is that a sig- 
nificant fraction (~ 30 %) of the metals ejected is retained 
in the galactic region when the moderate starburst SB1 is 
adopted. 



Figure 7. Map of the logarithm of the number density for model 
SB1 at four different times. Labels are in kpc. The gray-scale 
varies linearly from -5.42 (white) to 0.30 (black). Contours are as 
in Fig. 2 



two extreme models should bracket the reality. For this adi- 
abatic model the replenishing of the central region occurs at 
t ~ 85 Myr. 

Fig. 7c shows the ISM density at t = 100 Myr. The hot, 
rarefied bubble is almost totally shrunk; the hot gas mass is 
now only 4.4x 10 2 M Q (it was ~ 2.7x 10 4 M at t = 30 Myr). 
At the same time the adiabatic model contains ~ 2.6 x 10 4 
Mq of hot gas. The cold ISM continues to move ordinately 
toward the z-axis, and it encounters a weak accretion shock 
at R ~ 0.1 kpc. 

The density distribution at 200 Myr is shown in Fig. 
7d. No more hot gas is present (while in the adiabatic model 
Mhot ~ 4.8 x 10 3 Mq , and it is decreasing with time as the 
result of the numerical diffusion). The accretion shock has 
moved forward to R ~ 1 kpc, where the cold ISM is still 
accreting with v ~ 10 km s . The face-on surface density 
in the galactic region varies from 2.5 x 10 21 cm -2 at the very 
center, to 4 x 10 20 cm -2 at R = 2 kpc. The number den- 
sity in the central region is about 0.35 cm" 3 . As in model 
STD the central surface density is slowly increasing with 
time, approaching the critical value for the onset of effective 
star formation activity. Thus, also for this model, the sec- 
ular hydrodynamical evolution indicates the possibility of 
recurrent starburst episodes. The time between successive 



4.2 Model PEXT+SB2 

With this model we investigate the evolution of a galactic 
wind occurring in a galaxy immersed in a hot, tenuous ICM 
as described in section 2.1. All the other parameters are iden- 
tical to model STD. Fig. 8a shows the gas density at t — 30 
Myr. The superbubble has already blown out in the ICM, 
generating a complex filamentary structure. The fastest ma- 
terial penetrating in the ICM is moving with v w 2000 km 
s . Fig. 8b and 8c show the density at t = 60 Myr and at 
t = 200 Myr. The portion of the cold shell blowing out in the 
ICM is completely disrupted by the instabilities and spreads 
in a large volume, due to the high expansion velocities in the 
rarefied medium. At 200 Myr, the original ISM survives in a 
toroidal structure (1.5 < R < 8 kpc) on the equatorial plane. 
The inner edge of the cold gas is receding slowly (v ~ 20 km 
s _1 ) toward the center, while the outer portion is still ex- 
panding (v ~ 40 km s" 1 ). The cold gas starts to collapse 
toward the center, which is reached at t ~ 270 Myr, much 
later than in the previous models. The ISM column density 
increases more slowly than in model STD, and at t = 500 
Myr the central peak is only Eo ~ 2 x 10 20 cm -2 . Thus, in 
this case, the subsequent star formation episode might be 
delayed with respect to model STD. 

At 200 Myr the mass of gas present in the galactic region 
is MisM, g ai ~ 1.9 x 10 6 M Q , and about 1.5 % is hot (T ~ 10 6 
K). At the final time (500 Myr) we have Mi SM , g ai ~ 8.7 x 10 6 
Mq . The ejection efficiency is /ism = 0.31, much higher 
than in model STD because the absence of an extended en- 
velope of (relatively dense) cold gas. 

The mass of the metal-rich ejecta in the galaxy is 
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Figure 8. Map of the logarithm of the number density for model PEXT at three different times. Labels are in kpc. The gray-scale varies 
linearly from -7.27 (white) to 0.20 (black). Contours are as in Fig. 2 



M cjiga i ~ 5.8 x 10 3 M Q (1.8 x 10 4 M at 500 Myr) and 
fd = 0.83. These values are comparable to those found for 
model STD. However, we find that the hot gas has been 
severly contamined by the hot ICM, and Z < 0.05 for al- 
most all the hot ISM. The reason for this behaviour is the 
high temperature of the ICM, which greatly increases the 
importance of numerical diffusion. 

We estimate an upper limit for this effect, considering 
the first order upwind method (Roache 1972). The numer- 
ical diffusion coefficient is D upw ind ~ cA, where A is the 
zone size. The diffusion time is to = A 2 /D « 30 Myr (here 
A ~ 30 pc at R = z ~ 2.5 kpc and c is the ICM sound 
speed), so the numerical diffusion affects significantly this 
simulation, and this explains the very low values for Z. We 
conclude that for model PEXT we cannot calculate the en- 
richment process in a consistent way. For model STD, given 
the low temperature of the ISM (4.5 x 10 3 K), td is more 
than two order of magnitude longer, and the intrinsic diffu- 
sion is negligible. 

We note that the physical diffusion time scale, td = 
L 2 /D, where L is the typical length scale of the problem 
(L w 1 kpc), is very short: td ~ 10 2 — 10 3 yr. This is due 
to the high value for D « Ac, where A ~ 5 Mpc is the mean 
free path for the ICM (Spitzer 1962). However, even a small 
magnetic field reduces the mean free path to the order of the 
ion Larmor radius rr,. Only in this case we are allowed to 
consistently use the hydrodynamical equations. With A m tl 
the physical diffusion is effectively impeded. 

4.3 Model B+SB2 and BCOND+SB2 

Panel a of Fig. 9 shows the gas density of model B (sec- 
tion 2) at 30 Myr, just at the end of the starburst activity. 
The free expanding wind extends so far that almost all of the 
galaxy is devoid of the pristine gas. There is a radial gradient 
in the bubble temperature: along the 2-axis T ranges from 
~ 4 x 10 7 K close to the reverse shock to ~ 10 6 K behind the 



forward shock; a similar pattern is also present along the R- 
axis, although the temperatures behind the lateral shock are 
lower (~ 10 5 K) because of the lower velocity of the shock 
moving through the higher local ambient density. The aver- 
age density of the hot gas filling the bubble is ~ 10~ 4 cm -3 . 
The expansion velocity along the symmetry axis is ~ 300 km 
s _1 , and decreases toward the equatorial plane. The bubble 
accelerates as it expands through the decreasing ISM density 
profile, and the R-T unstable contact discontinuity generates 
relatively dense (n ~ 10~ 2 cm" 3 ) and cold (T ~ 10 4 K) fil- 
aments and blobs. Actually, denser structures can be seen 
on the z axis, but they are likely due to our assumption of 
vanishing radial velocity on the symmetry axis. In fact, cold 
gas deposited on this axis cannot be effectively removed, a 
well known shortcoming common to all 2D cylindrical simu- 
lations. Given the progressively increasing zones size with z 
and R, it is likely that the knots density is underestimated 
in our simulations, especially for the condensations far from 
the center. At t = 57 Myr (Fig. 9b) there is a large region 
essentially devoid of gas (n~5x 10~ 6 cm -3 ) surrounded 
by a very thick, low density shell (n ~ 10 -4 cm -3 ), with a 
temperature ~ 10 6 K. The external shock is rounder than in 
model STD, due to the lack of the collimating effect of the 
funnel along the z-axis (cfr. Fig. 2b). The dense and cold 
gas near the equatorial plane is already receding toward the 
center with a velocity ~ 30 km s" 1 , and a rarefaction wave 
is moving outward. The highest density in the shell is ~ 10 
cm" 3 on the equator, where the expansion velocity is 10 km 
s _1 . Apart the cold gas on the z axis, where the density 
reaches ~ 30 cm -3 , the densest filaments have n~3x 10~ 2 
cm -3 . At t ~ 75 Myr (not shown in Fig. 9) the inflowing 
cold gas reaches the center, where the density is still rather 
low (n ~ 10~ 2 cm -3 ), s -1 . After 106 Myr (panel c), the final 
time of this simulation, the hot gas is still expanding, but 
the central cold ISM is enterely collapsed, filling a region 
\z\ < 2.5 kpc, R < 6.5 kpc. The galaxy has thus recovered a 
cold ISM distribution similar to the original one. The ISM 
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Figure 9. Map of the logarithm of the number density for model B at three different times. Labels are in kpc. The gray-scale varies 
linearly from -6.31 (white) to 0.76 (black). 



Figure 10. Map of the logarithm of the number density for model BCOND at three different times. Labels are in kpc. The gray-scale 
varies linearly from -6.33 (white) to 0.95 (black). 



mass inside the galaxy is MisM, g ai ~ 10 Mq . Near the 
center the density reaches a few cm -3 . 

About the starburst ejecta, its largest content inside the 
galaxy is reached at t — 30 Myr, with M e j ;ga i = 3.6 x 10 
Mq . At this time < Z Ea x >= 2.4 x 10 -3 , and it decreases 
steadly with time. At t = 106 Myr < Z sa \ >= 1.1 x 10~ 3 
while the ejecta mass is Af c j, ga i ~ 10 5 Mq , about one tenth 
of the total gas lost by the massive stars. The ISM and 
metals ejection efficiency, estimated at t — 106 Myr, are 
/ism = 0.48 and / e j = 0.77 respectively. While the high / e j 
is consistent with the results obtained in section 3.3, the 
large value for /ism is striking when compared to model 



STD. However, this discrepancy reflects a deficiency in the 
definition of /ism, rather than a really different behaviour 
of the two models. As a matter of fact, model B recovers 
a 'normal' ISM before model STD! The discrepancy is due 
to several factors. First, in model B the massive dark halo 
is absent. The escape velocity is then quite low (a factor 
2-4 less than in model STD) and the gas residing at large 
radii becomes easily unbound. Second, the ISM distribution 
in model B is more peaked than in model STD, and the 
total amount of gas present in the numerical grid is lower 
(cfr. section 2.1). This, in turn, means that in model B the 
starburst provides more energy per unit gas mass than in 
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model STD. We believe that the difference in /ism between 
model STD and model B should be considered with some 
caution. In real galaxies, the gas at large radii, which is the 
source of the difference in /ism, can be removed by ram 
pressure and tidal stripping, processes not included in our 
simple models. Thus, the contribution of this gas to /ism is 
rather uncertain. 

The X-ray emission averaged < Z >x is much higher 
than < 2 ga i > and increases from < Z >x~ 0.06 at t = 30 
Myr up to < Z >x~ 0.2 at t > 50 Myr; at later times the 
bubble gas cools out of the X-ray temperatures and < Z >x 
drops to zero at t ~ 60 Myr. 

In Fig. 10 we show model BCOND, identical to model 
B but with the heat conduction activated. Again, panel a 
shows the density at 30 Myr. The superbubble is less ex- 
tended than in model B, because of the increased radiative 
losses in the conduction fronts. The temperature distribu- 
tion inside the superbubble is now rather flat near the equa- 
tor, but a negative gradient is present along the ^-direction, 
with the temperature in the range 10 6 < T < 10 7 K. As in 
the previous models, cold structures are present due to the 
R-T instabilities. The density of these structures is n ~ 3 
cm~ 3 , while the density of the hot gas is ~ 10~ 3 cm" 3 , one 
order of magnitude larger than in model B. This higher den- 
sity is due to the evaporation of the walls of the shocked ISM 
shell which 'feed' the inner region of the bubble. The expan- 
sion velocities are similar but lower than those of model B 
at the same time. Panel b shows the density at 56 Myr and 
can be compared with panel b of Fig. 9. The size of the su- 
perbubble remains smaller and the shape more elongated. 
The hot gas in the cavity is denser (n ~ 5 x 10 -5 cm -3 ) and 
slightly colder (T ~ 2.5 x 10 s K) than in the non-conductive 
case. The cold gas near the equatorial plane is receding to- 
ward the center, while the outer edge (where n ~ 10 cm~ 3 ) 
is still expanding. The densest filaments have n ~ 0.1 cm -3 . 
Panel c shows the gas flow at 105 Myr. The cold inflowing 
gas has just reached the center, much later than model B. 
In fact, the pressure drop of the hot gas is slower than in 
model B because the density of the hot gas is kept higher 
by the shell evaporation and by the slower expansion rate. 

At t = 125 Myr, the last time of this simulation, 
AfiSM.gai ~ 0.92 x 10 s Mq , not far from the initial value. 
However, only a fraction of the galactic volume contains a 
cold, dense ISM. In fact, roughly half of the galaxy is still 
filled with the rarefied gas of the cavity, now only moderately 
not ( T ~ !0 5 K) . At this time the ISM ejection efficiency is 
/ism = 0.32. 

The peculiar structure apparent on the simmetry axis 
(z > 20 kpc) at t — 105 Myr (Fig. 10c) is a numerical artifact 
depending on our treatment of the heat conduction. Colli- 
sionless shocks (for istance, in supernova remnants) do not 
show the hot precursor which would be expected (Zel'dovic 
& Raizer, 1966). This means that the plasma instabilities re- 
sponsible of the shock front formation also inhibit the heat 
flow through the front itself (Cowie 1977). To mimic this 
phenomenon in numerical simulations, the heat conduction 
coefficient must vanish at the shock front. To detect the 
shock front position on the computational grid is an easy 
task in ID simulations, but becomes rather cumbersome in 
two dimensions. Fortunately, the precursor length is rather 
short (shorter than the grid size) unless the upwind density 
is very low. We thus did not make any special treatment at 



the shock front. Effectively, the heat flux overruns the front 
only at late times, when the upwind density becomes rather 
low. However, this happens when the shock is well outside 
the galaxy, and our conclusions are not affected. 

The ejecta content inside the galaxy is Af c j = 1.7 x 10 5 
Mq after 30 Myr (< Z gal >~ 1.2 x 10~ 3 ) and decreases 
steadily down to Af cj = 0.77 x 10 5 M Q after 125 Myr, when 
< .Zgai >~ 8.3 x 10~ 4 . At t = 125 Myr we find / cj = 0.88. 

It is particularly interesting to investigate the X-ray 
emission for model BCOND, since now the numerical dif- 
fusivity does not affect the value of Lx and < Z >x (cfr. 
section 3.4). In fact, the thermal conduction naturally broad- 
ens the contact surfaces on lengthscales larger than thickness 
due to the numerical diffusion. The temporal variation of Lx 
is shown in Fig. 6. Lx is higher than in model B because 
of the emission arising in conduction fronts. The emission 
averaged abundance < Z >x ranges between 0.13 and 0.20 
for t < 70 Myr. As the energy input stops, the temperatures 
of the hot phase quickly drops and after t ~ 70 Myr no more 
X-ray emitting gas is present. The observable emission av- 
eraged temperature of the hot gas is < T 2 x 10 6 K 
for t < 30 Myr and drops quickly thereafter. 



5 DISCUSSION AND CONCLUSIONS 

The results presented here qualitatively confirm the conclu- 
sions of previous investigations on the effect of galactic winds 
in dwarf galaxies (e.g. MF). In general, it is found that the 
ISM is more robust than expected, and it is not disrupted 
even if the total energy input is much greater than the gas 
binding energy. In fact, the gas in the optical region of dwarf 
galaxies is only temporarily affected by the starburst, and 
the galaxy is able to recover a 'normal' ISM after a time of 
the order of 100 Myr from the starburst event, here assumed 
to be instantaneous. Our results agree well with the 'mod- 
erate form' of galactic wind dominated evolution of dwarf 
galaxies described by Skillman (1997). In Table 2 we sum- 
marize the values of the fraction of ISM and metal-rich ejecta 
that is lost by the galaxy. 

We find that the evolution of the ISM can be separated 
in two phases. The first one corresponds to the energy input 
period (which lasts 30 Myr in our models). During this phase 
the superbubble expands surrounded by a fragmented and 
filamentary shell of cold gas. The hot gas inside the bubble 
and the cold shell gas are in pressure equilibrium. The sec- 
ond phase starts when the energy input stops: the pressure 
of the hot bubble, still expanding along the polar direction, 
drops quickly. This causes the inner portion of the shell near 
the equator to collapse back toward the center, replenish- 
ing the galactic region with cold gas. The collapse is driven 
mainly by the pressure gradient, with the gravity being of 
secondary importance. 

The replenishment process occurs through inflow of cold 
gas moving parallel to the equatorial plane, thus resembling 
the inflows considered by Tenorio-Tagle & Munoz-Tunon 
(1997). However, the ram pressure associated to this flow 
in model STD, representative of all our models, is ~ 10~ 14 
dyn cm -2 , five orders of magnitude lower than those as- 
sumed by Tenorio-Tagle & Munoz-Tunon (1997). Evidently, 
if such massive inflows exist, they must have a different ori- 
gin. 
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Table 2. Global derived quantities. 
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-^ISM,gal 
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STD 
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0.46 
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PEXT 


0.31 


0.83 


1.9 x 10 6 


5.8 x 10 3 


3 x KT 3 


B 


0.48 


0.77 


9.3 x 10 7 


10 5 


1.1 x 10~ s 


BCOND 


0.32 


0.88 


9.2 x 10 7 


7.6 x 10 4 


8.3 x 10~ 4 



Quantities are calculated at t = 200 Myr, except for model B 
and BCOND, for which the final time is 106 Myr and 125 Myr, 
respectively, /ism is the ISM ejection efficiency defined in section 
3.2. / e j is the metals ejection efficiency defined in section 3.3. 
< -Zgal > is the average ejecta fraction in the galactic region (see 
section 3.3). Masses are given in Mq. 



We found that the central ISM reaches the critical col- 
umn density required for rapid star formation after 0.1-1 
Gyr from the starbust, the exact value depending on the 
galactic parameters, when a new starburst may start. This 
episodic star formation regime is necessary to account for 
the chemical evolution of BCD galaxies, and we have shown 
here that it is consistent with the hydrodynamical evolution 
of the ISM. 

Most of the metal-rich material shed by the massive 
stars resides in the hot phase of the ISM, and for powerful 
starbursts it is easily lost from the galaxy (Table 2). We 
estimate that a fraction of 0.5 — 0.9 of the total metal-rich gas 
is dispersed in the intergalactic medium when the starburst 
model SB2 is adopted. However, for moderate energy input 
rates (model SB1), only a small fraction (< 10 %) becomes 
formally unbound. In spite of the smallest / c j, model SB1 
has the lowest < Z sa i >, since the total amount of ejecta 
is a factor of 10 lower than the other models. Most of the 
ejecta material is pushed to large distance from the galaxy 
(several kpc), and its fate is uncertain, being subject to ram 
pressure and tidal stripping. These processes may effectively 
remove material loosely bound to the galaxy. 

There is some quantitative difference between our find- 
ings and those by MF. The generally lower / j found in our 
model is likely to be the result of our more extended gaseous 
halo. Their models, with a sharp truncation of the ISM, are 
similar to our model PEXT. The most striking disagree- 
ment is between model SB1, for which we find / c j = 0.003, 
and their model with M gas = 10 8 M Q and Li np = 10 39 erg 
s _1 which has / c j = 1. However, as explained in section 4.1, 
model SB1 suffers of some numerical extra-cooling; the same 
model without radiative losses gives / c j = 0.14, probably a 
more realistic value if thermal conduction is not effective. A 
comparison of our Fig. 7a-b with fig. 2a of MF (panel with 
the model M g = 10 s ; L^t, = 10, in their notations) dramati- 
cally shows the sensitivity of the superbubble dynamics (and 
size) on the ISM distribution. 

The cold gas replenishing the central region, from which 
a successive starburst may form, has been only slightly pol- 
luted by the massive star ejecta, with < Z >~ 4 x 10~ 4 — 2 x 
10~ 3 , or < Zq >~ 2 x 10~ 3 - 10" 2 Zo t@ , with the assump- 
tions described in section 3.3. Thus, many starburst episodes 
are necessary to build an average metallicity Z ~ 0.25 Zq, 
as determined for NGC 1569 from observations of emission 
lines gas (e.g. Kobulnicky & Skillman 1997). Many dwarf 
galaxies, however, are more metal poor (the abundance of 



IZwl8 is only 0.07 solar), and few bursts may be sufficient 
to produce the observed metallicity. 

The values listed in Table 2 demonstrate how the evo- 
lution of the ISM is not regulated by the ejection efficiency 
parameters /ism and / e j alone. For instance, model PEXT 
and BCOND have similar /ism and / e j, but very different 
ISM and ejecta masses. Conversely, the gas mass evolution of 
model SB1 and model B is comparable, despite of dissimilar 
ejection efficiencies. 

Comparing models STD and B, we found that a dark 
matter halo has little direct influence on the final behaviour 
of the ISM. However, dark matter and rotation determine 
the initial gas distribution, which is an important parame- 
ter for the flow evolution. For instance, the central region of 
model STD is refilled with cold gas more slowly than model 
B, a difference which reflects the different initial ISM distri- 
bution. 

In order to evaluate how the assumption of an instan- 
taneous burst influences our results, we run an additional 
model similar to STD, but with L inp = 1.128 x 10 40 erg 
s" 1 , M — 0.009 Mq yr" 1 . The energy and mass sources are 
now active for 100 Myr, so that the total energy and mass 
injected are the same as in STD. With this model, which in 
a simple way mimics the effect of a prolonged starburst, we 
wish to check the sensitivity of our results on the assumption 
of a instantaneous starburst adopted in models illustrated in 
the previous sections. We found that the general dynamics 
is similar to that of model STD and is not described here. 
The ejection efficiencies, again calculated at t = 200 Myr, 
are /ism = 0.055 and / e j = 0.40, almost identical to those 
found for STD (Table 2). Thus, it appears that the instan- 
taneous starburst hypothesis does not invalidate our general 
findings. 

We devoted section 3.4 to discuss the complex X-ray 
emission arising from starburst galaxies. We warn again that 
for model BCOND only we can calculate the X-ray quan- 
tities in a strictly consistent way, the other models having 
the contact surfaces numerically spread by the intrinsic dif- 
fusion of the hydrocode (although several physical processes 
are thought to produce similar effect, cfr. section 3.4). In 
model BCOND the thermal conduction generates physically 
broadened interfaces between hot and cold gas. We found 
that the X-ray luminosities in the ROSAT band (Fig. 7) 
are generally less than those estimated by Delia Ceca et al. 
(1996, 1997) for NGC 1569 and NGC 4449. It is thus sug- 
gested that a mass loading mechanism is at work in these 
systems (see also Suchkov et al. 1996). Moreover, the low 
abundances found in X-ray studies also indicate that ther- 
mal conduction (or some other process that mix cold and 
hot gas) is effective. 

Contrary to Suchkov et al. (1994) we found that the 
shocked ISM layer does not contribute appreciably to the 
X-ray luminosity for t < 30 Myr, a fact indicating that the 
origin of the X-ray radiation is model dependent. For Li np 
appropriate for typical dwarf galaxies we found that Lx is 
dominated by the ISM mixed with the shocked ejecta at the 
contact discontinuities. 
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